## This are customized functions for
## Robustness check
## comparing models including component measures
## Get robust se and pvalue
robust_se_pval <- function(dat, iv_total){
  G <- length(unique(dat$country_name))
  N <- nobs(iv_total)                     
  dfa <- (G/(G-1))* (N-1)/iv_total$df.residual
  clustered_vcov_ivtotal <- dfa*vcovHC(iv_total, type="HC0", cluster="group", adjust=TRUE)
  V.robust.se <- sqrt(diag(clustered_vcov_ivtotal))
  V.pvals <- coeftest(iv_total1, vcov= clustered_vcov_ivtotal)[,4]
  V.se.pval <- c(V.robust.se, V.pvals)
  return(V.se.pval)
}

## Get IV gof statistics
## texreg does not by default include 
ivreg_gof <- function(iv_total){
  tr <- texreg::extract(iv_total,include.rsquared=FALSE,include.adjrs=FALSE)
  iv_gof <- summary(iv_total,diagnostics=TRUE)$diagnostics[,4]
  names(iv_gof) <- paste(names(iv_gof), "test")
  if(length(iv_gof) <=3) {iv_gof = append(iv_gof, c("Weak instruments (I(1 + punish_percent_cumsum)) test"=NA,
                                                    "Weak instruments (e_h_polcon3_lag1) test"=NA), after=1)}
  iv_gof <- plyr::rename(iv_gof, replace=c("Weak instruments test" = "Weak instr. test: Demand",
                                 "Weak instruments (threatS) test" = "Weak instr. test: Demand",
                                 "Weak instruments (I(1 + punish_percent_cumsum)) test" = "Weak instr. test: Punished",
                                 "Weak instruments (e_h_polcon3_lag1) test" = "Weak instr. test: Elec. comp."),
               warn_missing=FALSE)
  tr@gof <- c(tr@gof, iv_gof)
  tr@gof.names <- c(tr@gof.names, names(iv_gof))
  tr@gof.decimal <- c(tr@gof.decimal, rep(TRUE, length(iv_gof)))
  return(tr)
}
